Structure and binding in crystals of cage-like molecules: 
hexamine and platonic hydrocarbons 
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In this paper, we show that first-principle calculations using a van der Waals density functional 
(vdW-DF), [Phys. Rev. Lett. 92, 246401 (2004)] permits determination of molecular crystal 
structure. We study the crystal structures of hexamine and the platonic hydrocarbons (cubane and 
dodecahedrane) . The calculated lattice parameters and cohesion energy agree well with experiments. 
Further, we examine the asymptotic accounts of the van der Waals forces by comparing full vdW-DF 
with asymptotic atom-based pair potentials extracted from vdW-DF. The character of the binding 
differ in the two cases, with vdW-DF giving a significant enhancement at intermediate and relevant 
binding separations. We analyze consequences of this result for methods such as DFT-D, and 
question DFT-D 's transferability over the full range of separations. 

PACS numbers: 61.50.Lt, 61.66.Hq, 71. 15. Mb 
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I. INTRODUCTION 

Understanding supramolecular structure and inter- 
actions is essential for understanding many biological 
processes^ Biological and other supramolecular com- 
plexes as polymers and overlayers are sparse matter, that 
is, they contain low electronic density in essential regions. 
The general lack of order in these systems prohibits pre- 
cise measurements of atomic structure and therefore chal- 
lenges development of theoretical methods. In turn, this 
makes transferable first principles schemes attractive; an 
accurate account of simple periodic structures permits 
accurate characterization to be made and reliable con- 
clusions to be drawn for more complex (non-periodic) 
systems. Molecular crystals and simple polymer crystals 2 
are ideal testing grounds for applications of first-principle 
descriptions of sparse matter. Unlike most sparse mat- 
ter, they constitute ordered systems and therefore leads 
to unambiguous comparison between theory and experi- 
ments. 

The dispersion forces underpin the cohesion of sparse 
matter. Modeling of sparse matter at the electronic level, 
therefore, requires that we take these effects into account. 
Traditional implementations of Density-Functional the- 
ory (DFT), which parameterizes the density functional 
within the local density approximation (LDA)^ or the 
generalized gradient approximation (GGA) 4 -^ have al- 
lowed for routine modeling of matter with dense electron 
distributions. However, the lack of non-local correlation 
and hence dispersion forces in these implementations, 
has effectively inhibited widespread use for general sparse 
matter systems. Computational methods able to predict 
molecular crystal structure and stability are also of value 
in the pharmaceutical industry, for example in mapping 
out competing crystalline phases for the fabrication of 
pharmaceuticals £ The electronic structure and response 
are often required to understand and compute properties 
of technological significance. Consequentially, a range of 
developments has aimed to extend density functional ap- 
proximations with an account of dispersive interactions 



and thus capability to address sparse matter challenges. 
Such approaches include non-local functional^— , inter- 
molecular perturbation theor y 13 : 14 , and the addition of 
a semi-emperical atom-based account of the dispersion 
forces (DFT-D) DFT-D has also been used to char- 
acterize molecular crystals . 19 ' 20 

We study molecular crystal using a van der Waals den- 
sity functional (vdW-DF) J£ The correlation part of this 
functional depends non-locally on the density and ac- 
counts for dispersion forces. vdW-DF has been used 
in a series of first-principle studies of sparse matter, 21 
yielding for example new insight to the twist of DNA, 22 
nanotube bundles^ 3 - water hexamersj 2 ^ metal-organic 
frameworks^ and the binding mechanisms at organic- 
metal interfaces 

In this paper, we first demonstrate that vdW-DF per- 
mits structural determination of molecular crystals and 
second asses the character of the attractive dispersion 
interactions. For the first, we calculate the lattice pa- 
rameters for crystals of the cage-like hexamine, dodeca- 
hedrane, and cubane molecules. For the second, we com- 
pare an atom-centered asymptotic- 1/r 6 approximation, 
frequently encountered in force-field methods ^ 30 : 31 and 
in DFT-D, with the correlation energy provided by vdW- 
DF. We also analyze consequences for DFT-D and ques- 
tion its performance at intermediate separations. The 
particular choice of molecules was motivated in part by 
their aesthetic appeal and in part for the property that 
their symmetric geometries permits a simple assignment 
and analysis of the strength of the asymptotic interac- 
tions even when expressed in terms of an atom-pair basis. 

This paper has the following plan: The next section 
presents the molecules hexamine, cubane, and dodeca- 
hedrane and their experimental crystal structures. The 
third section deals with the computational details of 
vdW-DF. The fourth gives the results for lattice param- 
eters and bulk modulus. In the fifth section, we compare 
the asymptotic l/r 6 -form of van der Waals interactions 
with the full non-local correlation for different molecules, 
and the vdW-DF potential energy curve with DFT-D cal- 
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FIG. 1: The molecular and crystal structure of hexamine, dodecahedrane, and cubane. The upper panels show the respective 
molecular structures. The mid panels give schematics of the crystal structures, where the molecular symmetry and orientation 
are highlighted by the use of a tetrahedron, dodecahedron, and a cube in place of the molecules. The lower panel shows, for 
hexamine and dodecahedrane, facets of the crystal structure and for cubane the orientation of the cubane molecule in the unit 
cell. Hexamine has a bcc unit cell, dodecahedron an fee unit cell, and cubane a rhombohedral unit cell, with equal angles a 
between the lattice vectors. 



culations for a cubane dimer. The final section holds our 
conclusions. 



II. PLATONIC MATTER 

Plato asserted that the basic building blocks of Na- 
ture were five geometrical structures, today known as 
platonic solids.™ These are the five convex polyhedra 
with all faces, edges, and angles congruent. On the 
molecular level, two of these solids have synthetic hydro- 
carbon (compounds of only carbon and hydrogen) ana- 
logues: cubane (CgHe), which corresponds to the cube, 
and dodecahedrane (C20H20), which corresponds to the 
twelve-faced dodecahedron. The fascination mathemati- 
cians held for these geometrical structures since antiquity 
echoed in the more recent struggle to synthesize their hy- 
drocarbon representatives. The highly strained bonds 
posed the main challenge. Obstacles were eventually 
overcome and in 1964 cubane was synthesized 33 followed 
by dodecahedrane in 1978<2i A third, tetrahedrane exists 



only within a larger chemical structure 3 ^ and no crys- 
tallographic characterization exists. Tetrahedrane does 
therefore not represent a good testing ground for vdW- 
DF and we instead study the molecular crystal of hexam- 
ine (C6H12N4). Although not a hydrocarbon, it shares 
several features with the platonic hydrocarbons: it is or- 
ganic, its non-hydrogen atoms form a cage, and this cage 
has the symmetry of platonic solid tetrahedron. The top 
row of Fig. [T]shows the molecular structures (from left to 
right) of hexamine, cubane and dodecahedrane, ordered 
according to the ascendancy of their platonic analogues. 

The second and third row of Fig. Q] show the crystal 
structures. The first column shows hexamine, the mid, 
cubane, and the final, dodecahedrane. The crystal struc- 
ture of hexamine forms a bcc with a single molecule in 
each unit cell and I43m symmetry. Hexamine has been 
used as a model system in numerous studies 3 ^— and the 
crystal structure was determined as early as in 1923. 39 
The crystal structure of cubane is rhombohedral with 
equal external angles and a single molecule in each cell. 
The molecule is oriented according to the R3 space group; 
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a rotation about the [111] axis specifies the configuration 
of the molecule* 3 -! The crystal structure of dodecahedrane 
is fee, also with a single molecule in each cell and Fm3- 
space group 42 

Hexamine finds numerous industrial uses 4! For in- 
stance, it serves as a component in fuel tablets and as 
an antibiotic4 2 - The platonic hydrocarbons find only hy- 
pothetical applications; cubane has been identified as a 
potential high-energy fuel and explosive. 43 

In addition to being well-studied crystals, especially 
hexamine and cubane, these crystals make attractive 
testing grounds for vdW-DF for two more reasons. First, 
their simple structures allow for brute-force determina- 
tion of lattice parameters. Within the crystal symmetry, 
this determination corresponds to mapping out the po- 
tential energy of a single parameter for hexamine and do- 
decahedrane, and three for cubane. The brute-force ap- 
proach facilitates post-processing analysis, such as com- 
putation of bulk modulus^ 4 - and makes it easier to eval- 
uate how the choice of exchange in DFT affects crystal 
structure and cohesion energy. Second, the high symme- 
try of the molecules reduces the large set of atom-to-atom 
C*6 coefficients to only a few equivalent values. This re- 
duction simplifies the comparison of full vdW-DF calcu- 
lations with atom-based asymptotic accounts of the van 
der Waals forces. 



III. COMPUTATIONAL DETAILS 

The crystal structure and bulk modulus of hexamine, 
cubane, and dodecahedrane are determined with DFT 
using a nonlocal density functional called vdW-DF. The 
traditional semi-local GGA for exchange-correlation pro- 
vides accurate bond lengths and charge density n(r), 
but fails to capture correlated motion of separated elec- 
trons: the long-range dispersion forces. vdW-DF in- 
cludes these correlations and can therefore account for 
the structure and cohesion of sparse matter. Since de- 
tails of the functional and its implementation are given 
elsewhere ) 10 ! 21 ' 45- — we focus mostly on computational 
steps specific for determination of crystal structures. 

The non-local correlation of vdW-DF takes the form 
of a double-space integral, 

Ef[n] = - f dr [ dr'n(r) <p(r, r')n(r') , (1) 
2 Jvo Jv 

over an interaction kernel </)(r,r'). Here Vb denotes the 
central unit cell and V (formally) the entire space. The 
kernel can be tabulated in terms of two parameters d and 
d' related to the local response qo(r) and spatial separa- 
tion |r— r'| by d = go (i") |r — r'| and d' = <7o( r ')l r— r '|- The 
remaining part of the exchange-correlation functional of 
vdW-DF consists of the exchange part of revPBE 4 ^ and 
the correlation of LDA: 

£VdWDF = E LBA + ^revPBE + £ nl [n] _ (2) 



The total energy functional of vdW-DF, E vdWDF [n], 
also includes the standard elecrostatic and kinetic-energy 
terms within the Kohn-Sham scheme.— It is convenient 
to write this energy in terms of a semi- local part -BoM 
containing all but the non-local correlation, so that 
E vdWDF [n] = Eo[n] + Ef[n]. For input charge density 
n(r), we use the result of semi- local calculations with 
the PBE 5 flavor of GGA. We will refer to calculations 
with the PBE flavor of GGA as DFT-GGA. The charge 
density could also have been obtained within fully self- 
consistent vdW-DF4£i££ However, the two-step non-self 
consistent procedure introduces only a slight approxima- 
tion. Previous studies have documented that for systems 
with small charge transfer, binding energies of non-self 
consistent vdW-DF only differ by tiny amounts from fully 
self-consistent energiesj 45 ' 51 

To speed up evaluation of the non-local correlation, we 
introduce a radius cutoff based on the decay of van der 
Waals forces at large separations. With this cutoff, the 
kernel takes form 

Ef[n]&- ( drn(r) [ dr' <£(r, r')n(r') . (3) 

2 Jy J\t-t'\<R 

In the above expression, we see that the CPU-cost 
for evaluating the non-local part of vdW-DF goes as 
R 3 0(Vo)- Thus, for large or periodic systems the compu- 
tational costs increase linearly with system size. To cut 
computational costs further, we introduce an extra radius 
cutoff corresponding to the separation between dense and 
sparse sampling of the charge-density grid. We note that 
Ref. |50| reports a more elaborate scheme which consid- 
erably reduces CPU-costs, yet our simple measure was 
sufficient for our non-self consistent calculations as the 
underlying DFT-GGA calculations of electronic density 
dominated time consumption. 

The use of revPBE for exchange in vdW-DF was mo- 
tivated by the fact that this exchange functional ex- 
cludes unphysical binding effects at large distances ! 46 ' 52 
For a range of systems, vdW-DF overestimates binding 
separations! 21 ' 46 Several studies indicate that this dis- 
crepancy can be attributed to the details of the exchange 
functional ! 21 ' 51 ' 53 " — Puzder et al 5 - 3 - have demonstrated 
that replacing Hartree-Fock exchange with revPBE im- 
proves binding separations for benzene dimers. Gulans 
et have shown that for a selected range of molecular 
complexes the PBE exchange functional improves bind- 
ing energies. We furthermore illustrate the sensitivity to 
exchange by including results based on use of an alterna- 
tive vdW-DF(PBE), where revPBE exchange has been 
replaced by that of PBE. We do not argue for replacing 
vdW-DF with vdW-DF(PBE), instead we simply explore 
consequences of a different account of exchange 4 s - 

To calculate crystal parameters and cohesion energy, 
we minimize the potential energy. In many respects this 
vdW-DF structure determination is similar to those in 
Refs. HUHH The potential energy is given by the 
difference between the total energy of the full crystal and 
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a reference energy for a system of isolated molecules 

E coh {a, {a0}) = E vdWDF (a, {a9})-E vdWDF {a -> oo, {afl}) 

(4) 

In the above equation the curly brackets are specific for 
the cubane crystal as it depends on three rather than one 
parameter. In the reference calculation corresponding to 
the reference energy, i5 vdWDF (a — > oo, {a.9}), the semi- 
local part is obtained somewhat differently from the part 
containing the non-local correlation. For the former, we 
effectively isolate the molecules in our periodic boundary 
calculations, by using a unit cell of doubled size in all di- 
rections. This measure secures negligible charge overlap 
between the molecules in the supercells. For the latter 
(nonlocal) part we restrict the integral of Eq. [3] to the 
central supercell to avoid coupling between the enlarged 
unit cells. Hence, only non-local correlations within the 
molecule contribute to the reference energy. 

To enhance accuracy in the evaluation of the non-local 
part of the potential energy, we systematically cancel a 
small, but noticeable, grid dependence in the evaluation 
of Eq. [3] This cancellation is performed by making sure 
to use the same FFT grid spacing in the reference calcu- 
lation as in the main calculation. Furthermore, we make 
sure to place the isolated molecules in the same relative 
configuration to the underlying grid in the reference cal- 
culation as in the main calculation. We thus perform an 
additional reference calculation for every molecular con- 
figuration investigated. These measures have been used 
to secure a high accuracy of the non-local part of vdW- 
DF in several earlier studies.— i 46 i 47 i 56 

We map the potential energy landscape by varying the 
lattice parameters of the molecular crystals within the 
experimental crystal symmetry. The stiff cage-molecules 
allow us to keep the internal coordinates of the molecules 
frozen for all configurations. The molecular structures 
arc determined in isolation using the PBE flavor of GGA. 
The resulting bond lengths will be compared with ex- 
perimental data in the next section to verify the utility 
of conventional DFT-GGA for the internal structure of 
strained molecules. 



TABLE I: Experimental and calculated bond lengths of hex- 
amine, cubane, and dodecahedrane. The calculations were 
done with the PBE flavor of GGA. I denotes the C-C bond 
length for cubane and dodecahedrane and the C-N bond 
length for hexamine. Iqh denotes the carbon-hydrogen bond 



length. 


Parameter 


Hexamine 


Cubane 


Dodecahedrane 


l[k] 

r xp [A] 

*CH [A] 


1.472 
1.476" 

1.101 
1.088" 


1.566 
1.562 1 ' 

1.095 
1.097 6 


1.549 
1.544 c 
1.100 


a Ref. 
6 Rcf. 
c Ref[ 


58. 



TABLE II: The vdW-DF prediction of lattice parameters, co- 
hesion energy and bulk modulus for the crystals of hexam- 
ine, cubane and dodecahedrane compared with an alternative 
vdW-DF(PBE) based on PBE-exchange and with experimen- 
tal values. The experimental lattice parameters are based on 
low temperature measurements, except for dodecahedrane. 





Parameter 


vdW-DF 


vdW-DF (PBE) 


Exp. 


Hex. 












a (A) 


7.14 


6.93 


6.910" 




E coh (eV) 


-1.0127 


-1.427 


-0.827 6 




Bo(GPa) 


10.0 


14.0 


7.0 C 


Cub. 












a (A) 


5.45 


5.25 


5.20 d 




a 


73 


72.5 


72. 7 d 







47.5 


46.5 


46 d 




E coh {eV) 


-0.77 


-1.15 


-0.857 c 




Bo(GPa) 


7.2 


14.8 




Dod. 












a (A) 


10.92 


10.56 


10.60 / 



E coh (eV) 
B {GPa) 



-1.46 
12.2 



-2.06 
18.6 



a Ref. 58. 6 Rcf. 63- c Ref.|62. d Ref.[5|. e Ref.|63|. ^Ref. 64. 



The electronic-structure calculations rely on the plane- 
wave code DACAPD 57 using ultra-soft pseudopotentials. 
In combination with a separate reference calculation as 
previously discussed, we secure the convergence of the 
non-local correlation by specifying an FFT-grid spacing 
less than 0.13 A. This spacing leads to an effective plane- 
wave energy cutoff of at least 500 eV. For all crystals, we 
set the Monkhorst-Pack k-sampling to 4 x 4 x 4. 



IV. STRUCTURE DETERMINATION 
A. Molecular structure 

Table Q] shows the calculated molecular structures. 
It demonstrates DFT-GGA can account for the in- 
tramolecular bonding even in the highly strained cubane 
molecules. The resulting bond lengths differ by less than 
1 % from the experimental values. In first principle stud- 
ies of molecular crystals, accurate determination of lat- 
tice parameters require accurate account of constituent 
molecules. 

Unlike many calculations, where determination of crys- 
tal structure often starts from the structure of the indi- 
vidual molecules, most experiments resolve the molecu- 
lar structure by looking at the diffraction pattern of a 
full molecular crystal. For our purposes of testing the 
first-principle vdW-DF method, this is fortunate, as ef- 
forts to characterize molecular structures also generate 
an abundance of experimental data on molecular crystal 
structures. 
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FIG. 2: The potential energy curves for the crystal of hexamine and dodecahedrane (left panel) and corresponding contour plots 
for the crystal of cubane (right panel). The left panel shows potential energy where the curves are normalized separately so that 
the experimental lattice parameter and the calculated cohesion energy equals unity. The solid curve represents dodecahedrane 
and the dashed hexamine. The difference in the optimal a/o eX p value can be attributed to the experimental lattice parameter 
being measured at low temperature for hexamine, but not for dodecahedrane. The right panel shows spline-interpolated contour 
plots of the two-dimensional intersection of the three-dimensional potential energy landscape. The main figure corresponds 
to the optimal value of the internal angle 9, while the insert corresponds to the optimal value of the unit cell length a. The 
pronounced asymmetry of the curves and contours reflect the hard wall provided by Pauli repulsion. 



B. Crystal structure 

Figure [2] shows the binding curves and contours found 
by varying the molecular crystals' lattice parameters 
identified in Fig. Q] The simple crystal structures of 
hexamine and dodecahedrane give a one-dimensional po- 
tential energy landscape, while the cubane crystal has a 
three dimensional one. In the right panel of Fig. [21 we 
display the aa and the a6 intersections, for the optimal 
values of 9 and a. The curves exhibit a pronounced asym- 
metry around their minimum. This asymmetry arise be- 
cause the van-der Waals attraction is much softer than 
the kinetic-energy repulsion. 

Table. |TT] contains the calculated results and experi- 
ment values. Standard vdW-DF performs well both for 
lattice parameters and cohesion energy. If not directly 
available, we obtain experimental cohesion energies by 
correcting for gas phase and vibrational contributions to 
the enthalpy of sublimation, using the method described 
in Refs. EH 

The bulk modulus is obtained with use of polynomial 
interpolation according to the scheme of Ziambaras and 
Schroder. 4 "* The required polynomial were constructed 
using data from selected one-dimensional deformations. 67 
For hexamine, where the experimental bulk modulus is 
available, the computed value show fair agreement with 
the experimental value. There is also a trend for bigger 
cage molecules to have a larger bulk modulus. We at- 
tribute this trend to the fact that for bigger molecules a 
smaller relative part of the unit cell consists of soft intra- 
molecular regions. Therefore, as the relative unit cell di- 
mensions change, the distance between (stiff) molecules 



changes more for big molecules than for small molecules. 

For all three crystals, we find unit-cell volumes some- 
what larger than the experimental ones. Similar over- 
estimations have also been encountered in previous 
studies ! 21 ! 46 ' 55 The alternative choice of PBE as exchange 
functional, vdW-DF(PBE), influences results substan- 
tially. On one hand, it improves lattice parameters, al- 
most to level of standard DFT for intra-molecular bonds. 
On the other hand, the value of cohesive energy and bulk 
modulus worsens. These results signal that the main dis- 
crepancy between experiments and vdW-DF stems from 
the specific form of semi-local exchange! 21 ' 5 3 ~ — 

We note that an ab initio study of cubane has previ- 
ously been performed at the LDA levels Being a local 
functional, LDA has no physical basis for the van der 
Waals binding that provides the cohesion of this molec- 
ular crystal. The spurious LDA binding arises from an 
unphysical accounts of exchange! 46 ' 52 Once the molecu- 
lar crystals are investigated with DFT-GGA, which has 
an improved account of exchange, the binding essentially 
vanishes^ 



V. ASYMPTOTIC PAIR POTENTIALS VERSUS 
NON-LOCAL CORRELATION 

The asymptotic van der Waals interactions between 
two atoms or molecules goes as the simple power law 
Cg/r 6 , where Cq gives the strength of the interaction. 
This familiar result can be derived from second-order per- 
turbation theory, 69 or from an analysis of the shifts in the 
zero-point motions of the electron^ A common strategy 
in force-field methods and empirical extensions of DFT 
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d-d [k] 



than to the surface atoms. Within vdW-DF, Kleis et a@£ 
demonstrated that for interactions in nanotube bundles, 
the force stems primarily from the electron tail around 
the nanotube, and that as the tubes get closer, higher or- 
der moments dominate over the asymptotic interaction. 
Part of this enhancement can be interpreted as an image- 
plane effect. 

Here we investigate whether an atom-centered 1/r 6 
form is a good approximation for the non-local corre- 
lation of molecular dimers. As our argument is based 
on the asymptotic vdW-DF account of van der Waals 
forces, we also discuss for these molecules the accuracy 
of the asymptotic account. In the first subsection, we 
compare the non-local correlation of vdW-DF for dimers 
of hexamine, cubane, and dodecahedrane with its corre- 
sponding atom-centered pair potentials. We document 
a significant enhancement of the non-local correlation at 
short (binding) separations and at intermediate separa- 
tions, one to three Angstrom beyond typical binding sep- 
arations. In the second subsection, we discuss the accu- 
racy of our Cq coefficients. In the third, we investigate 
consequences of this result, in particular for the use of 
DFT-D. We will argue that although standard DFT-D 
methods can provide good descriptions of both short and 
asymptotic separations, their asymptotic atom-centered 
form does not describe the enhancement of correlation at 
intermediate separations exhibited by vdW-DF. 



FIG. 3: Comparison between the non-local correlation of 
vdW-DF and atomic pairs potentials generated with asymp- 
totic vdW-DF for a dimer of hexamine, dodecahedrane, and 
cubane (from top to bottom) in configurations correspond- 
ing to nearest neighbors in their respective crystal. The 
dashed curve gives the APP-vdWDF result, while the solid 
curve gives non- local correlation of vdW-DF: AEf(d) = 
Ef(d) - Ef(d -¥ oo). The horizontal axis gives the differ- 
ence from the vdW-DF crystal binding separation do, which 
respectively takes the values 6.18, 5.45 and 7.72 A. The differ- 
ence between the two curves demonstrates the enhancement 
of non-local correlations over the vdW-DF asymptotic atom- 
based account at relevant binding separations and intermedi- 
ate separations. 



is to adopt such an asymptotic form at all separations in 
terms of atom-centered pair-potentials (APP). However, 
we can not take for granted that the asymptotic behavior 
should hold for separations closer to that of intramolec- 
ular binding. On the contrary, because van der Waals 
forces arise from correlated motion of electrons and not 
from the atomic nuclei, several mechanisms affect and 
enhance the interaction at short and intermediate sepa- 
rations: higher order moments contribute, polarizability 
changes as charge is distorted and finite-k dispersion of 
the electronic response becomes important. Zarcmba and 
Kohn 71 considered adsorption of noble atoms on surfaces 
and documented a significant enhancement of dispersion 
energies over an atom-centered account; their asymptotic 
1/d 3 form use the distance to an image plane d, rather 



A. Asymptotic vdW-DF 

For the the asymptotic van der Waals forces, the Cq 
coefficient between two fragments, A and B, can be com- 
puted from the general formula 



C, 



AB 



du aA(iu)aB{iu) , 



(5) 



where a(w) is the polarizability of the fragment. In order 
to calculate the Cq coefficients we approximate a(u>) with 
the local external-field susceptibility of vdW-DF, 



vdW-DF 



ru(r) 



[9 go (r) 2 /8^] 
and a polarizability given by 



(6) 



, vdW-DF 



(w) = J d 3 r 



3 "X v / W -° F (^,r). (7) 



To generate vdW-DF based atom-centered pair poten- 
tials (APP-vdWDF), we first partition the full charge 
density n(r) = J2i n i( r ) among the atoms of the 
molecules with aid of Bader analysis.— >S Based on this 
charge partition, we calculate the atom-to-atom Cq coef- 
ficients using Eq. ([5]). Initially, this procedure generates 
N 2 different C§ coefficients for a molecular dimer of N 
atoms per molecule. For a dodecahedrane dimer, we get 
as much as 40 x 40 = 1600 coefficients. Fortunately, be- 
cause of the high symmetry of the isolated molecules, 
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this number reduces to only three equivalent values for 
cubane and dodecahedrane: Cq~ c , C,?~ H , and Cf~ H . 
For hexamine, the extra nitrogen atoms lead to six co- 
efficients. As noise in the electronic density affects the 
value of the coefficients, we average over a large set of 
equivalent values to obtain the final values. 



TABLE III: Computed values of the C& coefficients (Hartree 
atomic units) for different pairs of atoms within respective 
molecules, using the asympotic form of vdW-DF with charge 
density as partitioned with a Bader analysis. The figure also 
shows the molecule-molecule Ce coefficients obtained with 
the Andersson-Langreth-Lundqvist (ALL) scheme,— and the 
molecule-molecule Ce coefficents obtained with the use of pa- 
rameters given in DFT-D schemes.— 



yivdWDF 


hex 


cub 


dod 


C-C 


4.44 


13.5 


11.2 


C-H 


4.04 


6.87 


6.45 


H-H 


3.98 


3.72 


4.07 


N-N 


32.6 






N-C 


12.0 






N-H 


11.0 






mol-mol (vdW-DF) 


3470 


1990 


11300 


mol-mol (ALL) 


3270 


1940 


10200 


mol-mol (Wu) 


4340 


2600 


16300 


mol-mol (Grimmo) 


4000 


2630 


16400 


mol-mol (Jurecka) 


4790 


3120 


19500 



The upper part of Table Hill shows the C§ coefficients 
calculated within asymptotic vdW-DF, both as parti- 
tioned according to the Bader analysis and as evaluated 
for the entire molecule. The calculated coefficients per 
atomic pair deviates much from a naive assignment of 
the full coefficient of the molecule according to the num- 
ber of valence electrons of the underlying atom. In such a 
scheme the C-C coefficient would be 16 times larger than 
the H-H coefficient, while in fact, for cubane and dodec- 
ahedrane, it is merely three-four times stronger. This re- 
sult can be attributed to the relative stronger response of 
the low-density regions surrounding the hydrogen atoms 
asgo(r) <x [n(r)] 1 / 3 (in the homogeneous limit), and these 
areas dominate Efr. The somewhat anomalous values for 
hexamine can be attributed to our Bader analysis scheme 
partitioning a significant portion of the charge density 
near the carbon atoms to the centrally located nitrogen 
atoms. 75 For cubane and dodecahedrane the partitioning 
was similar, and the differences in atomistic Cg values 
show that they are influenced by their local enviroment. 

Having generating C6 coefficients appropriate for a 
comparison between the asymptotic account and the 
full correlation of vdW-DF, we study dimers of hexam- 
ine, cubane and dodecahedrane at different separations. 
We choose orientations that are given by the nearest- 
neighbor configurations in the crystals. The total asymp- 



totic non-local correlation of APP-vdWDF reads 

K PP -2^2^ | r -r'l ' ( ' 

i j % V 

where i and j label atoms in separate molecules of the 
dimer. 

Figure [3] shows results for the full non-local correla- 
tion of vdW-DF (solid curve) and that of APP-vdWDF 
(dashed curve). At large separations these two curves 
converge. In contrast, they differ significantly at relevant 
binding (short) separations and at intermediate separa- 
tions. At these separations the non-local correlation is 
almost double as large as that of APP-vdWDF, for both 
hexamine and cubane, while for the biggest molecule, do- 
decahedrane, the difference is smaller. For all dimers, we 
also need to go to relatively large separations to recover 
asymptotic values for the non-local correlation. 

As vdW-DF is an ab initio functional, based on a set of 
exact sum rules^ the strong enhancement of non-local 
correlations at short and intermediate separations indi- 
cates that the asymptotic form neglects important contri- 
butions. It also highlights that in constructing modified 
semi-empirical van der Waals functionals) 11 ' 12 a fitting 
of the asymptotic functional to Cg coefficients does not 
guarantee an accurate description at short binding sepa- 
rations. 



B. Comparison of C'e coefficients 

The lower part of of Table IIII1 gives the molecu- 
lar Cq coefficients calculated with asymptotic vdW-DF, 
the Andersson-Langreth-Lundqvist (ALL) scheme^ and 
computed using the atomic coefficients of Wu el al^ 
Grimme, 17 and Jurecka et alJ£ The ALL and vdW- 
DF give coefficients smaller than that used in DFT-D 
schemes. 

We expect that Grimme and Wu provide good values 
for molecular Cq coefficients, because their underlying 
atomistic coefficients were fitted to reproduce a range 
of accurate molecular C§ coefficients calculated from ex- 
perimental molecular polarizabilities (Ref. [l6| and refer- 
ences therein). The coefficients of Jurecka^ give some- 
what larger molecular C§ coefficients. This comes from 
the use of Slater-Kirkwood average 7 ^ for Cq coefficients 
between different atomic species, while keeping those of 
Grimme for indentical atomi species (the Cq coefficients 
of Grimme are optimized for a different average. 

Asymptotic vdW-DF and ALL likely underestimates 
the Cq coefficients for these molecules; they are 20-30 % 
smaller than that used in DFT-D methods. For the sim- 
ilar ALL scheme, Ref. [7?] reports an underestimation of 
Cq coefficients for larger molecules, in particular for ben- 
zene and C60. A difference between the ALL scheme and 
asymptotic vdW-DF is that for the former a hard cutoff 
accounts for plasmon damping, while for the latter, the 
local response qo(r) provides a smooth cutoff. There is 
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good consistency between the two methods. Both meth- 
ods also assume a local, scalar, relationship between the 
applied and the full electric field, which is an approxima- 
tion for finite-sized objects>2i We speculate that this ap- 
proximation contributes to the underestimation of Cq co- 
efficient for the investigated, relatively large, molecules. 



TABLE IV: Atomistic Ce coefficients for cubane used in APP- 
vdWDF (calculated with asymptotic vdW-DF and charge 
density partitioned with Bader analysis), and coefficients used 
in DFT-D methods. 



r.vdWDP 


vdW-DF 


Wu° 


Grimme 


Jurecka c 


C-C 


13.5 


22.06 


28.3 


28.3 


C-H 


6.8 


7.89 


5.01 


8.82 


H-H 


3.7 


2.83 


2.75 


2.75 



°Ref. QJ 
''Rcf. 
c Ref. 



Tabic HVl shows atomistic C§ coefficients for cubane as 
calculated with asymptotic vdW-DF, and given by Wu, 15 
Grimme*^ and Jurecka^ for use in DFT-D schemes. 
vdW-DF weights the relative response of the carbon less 
than that of the hydrogen, compared with the values 
of DFT-D i£i This property could relate to the above 
mentioned approximate treatment of electrodynamics. It 
could also relate to the carbon atoms being located some- 
what inside the molecule, having a different local charge 
density and responding less to external fields than an 
atom on the exterior would; in contrast, DFT-D does 
not discriminate between atoms at different locations. 



C. Consequences for atom-based pair potentials: 
the binding curve of cubane 

The vdW-DF results presented in the first subsection 
shows that a simple asymptotic account only partially 
captures correlation effects at short and intermediate sep- 
arations. As DFT-D use such an asymptotic form to 
describe non-local correlations, this result stands in ap- 
parent contrast to the many successful applications of 
DFT D 16 - 18 - 78 - 80 

To understand consequences of our result for meth- 
ods such as DFT-D, we must first consider other effects 
that could contribute to the difference between APP- 
vdWDF and the non-local correlation of vdW-DF. In 
vdW-DF, correlations are described by ££ DA + Ef and 
hence E™ 1 also accounts for semi-local correlations ! 10 i 46 
Second, we must consider the specific designs of actual 
DFT-D schemes, because these could counteract the lack 
of enhancement of non-local correlations. To describe 
exchange-correlation, DFT-D combines the asymptotic 
atom-centered form with a semi-local GGA account. It 
also introduces fitting parameters to be used in combi- 
nation with a specific GGA flavors. 



The top panel of Fig. [4] shows that gradient corrections 
does not account for the difference between the full vdW- 
DF and the APP-vdWDF results. For APP-vdWDF, 
we can combine the purely semi-local correlation of PBE 
with APP-vdWDF (in a new description APP-mod), 



E 



vdWDF 
APP-mod 



M = < P r F (rf) + A£* BE (d) - A£ c LDA (d) , 

(9) 

to assess the magnitude of purely semi-local corrections 
relative to the difference between vdW-DF and APP- 
vdWDF. In this APP account, we have introduced the 
LDA and PBE terms: AE c (d) = E c (d) - E c (d -> oo). 

We focus our discussion on the cubane dimer. The 
lower thin solid curve gives APP-mod; the thick solid 
curve gives the non-local correlation of vdW-DF. The 
curves are shown as a function of the intermolecular dis- 
tance d, with do indicating the vdW-DF binding distance 
in the crystal. At short separations, the thin curve lies 
closer to the thick curve than the corresponding APP- 
vdWDF result (thick dashed curve) . Thus, some of the 
difference between APP-vdWDF and the non-local cor- 
relation arises from a lack of semi-local correlation con- 
tributions in APP-vdWDF. 10 However, even with this 
inclusion the difference is still significant, and at inter- 
mediate separations it remains undiminished. 

The middle panel of Fig. @] details the effects of using 
the semi-empirical fitting of DFT-D and the larger Cq 
coefficients. We compare the binding curve for cubane 
obtained by use of vdW-DF with the binding curve ob- 
tained with DFT-D calculations. We select the schemes 
of Grimme^ and Jureckaj^ as these provide parame- 
ters for the PBE flavor of exchange-correlation, which is 
available to us^ 

The DFT-D scheme of Grimme scales the strength 
of the dispersive interaction to a particular semi-local 
exchange-correlation to achieve good performance at 
short separations. For the PBE flavor of GGA, the cohe- 
sion energy of a dimer reads 



^DFT-D.G _ pPBE 
^coh — ^coh " 



PBE \ " rij , 

b 6 / J J PBE.G' 



r*3 
°6 



(10) 

where sq — 0.7. This scheme therefore sacrifices the 
asymptotic description in favor of the a good description 
of binding separations. The scheme of Jurecka instead 
adjust parameters of the damping function to the flavor 
of semi-local exchange-correlation, 



E 



DFT-D, J 
coh 



pPBE 
^coh 



E 



/pBE,j(l r i- r jl)l 



(11) 

and ensure that /p BE j — > 1 for large separations. The 
DFT-D scheme can therefore, in principle, describe both 
asymptotic and binding separations. However, the ques- 
tion remains on how it performs for systems where char- 
acteristic separations lies between these two limits. 

The binding curves of vdW-DF(revPBE) (the solid 
curve) and for the DFT-D scheme of Grimme (PBE) (up- 
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d - do[A] 



FIG. 4: Comparison between different accounts of non-local 
correlation for a dimer of cubane. In the upper panel, the 
dashed curve gives APP-vdWDF. The thin solid curve gives 
the sum of APP-vdWDF and the correlation of PBE. The 
thick solid curve gives the non-local correlation of vdW-DF. 
In the middle panel, the thick curve gives the cohesion en- 
ergy using vdW-DF, £'^ VDF (d), the upper (lower) dash- 
dotted curve gives the DFT-D binding curve as given by 
Grimme, E° l T -°' G (d) (Jurecka, £ c D h T " D,J (d)). The lower 
panel gives the corresponding non-local correlation of vdW- 
DF and Grimme (Jurecka). The thin solid curve gives the dif- 
ference between exchange of revPBE and PBE. The horizon- 
tal axis gives the difference from the vdW-DF crystal binding 
separation do for cubane. 



per dash-dotted curve) and Jurecka(PBE) (lower dash- 
dotted curve) indicates an underestimation of DFT-D 
at intermediate separations. The minimum at negative 
d — do shows that use of DFT-D would improve lattice 
constants over vdW-DFi^ The two DFT-D schemes give 
quite differing binding energies. Both results show that 
a binding energy at the same magnitude as vdW-DF can 
be achieved, even with an atom-based asymptotic form 
of the attractive potential. The energy of vdW-DF is sig- 
nificantly larger at intermediate separations than that of 
Jurecka despite that Cg coefficients of asymptotic vdW- 
DF are underestimated (while those of Jurecka are likely 
to be somewhat overestimated) and despite that the cor- 



responding DFT-D scheme binds stronger than vdW-DF. 

The lower panel of Fig. [4] shows that the difference at 
intermediate separations comes primarily from the vary- 
ing accounts of correlation. The solid curve shows the 
non-local correlation of vdW-DF. The dash-dotted up- 
per (lower) curve shows 

£DFT-D,G(J) = fij fir- — r-H ^ 

^app — /—I • / PBE,G(J)vl l1 jl '| r . -r-| 6 

ij 3 

+ AE PBE {d) - A£ c LDA (d) , (12) 

which for DFT-D corresponds best to the non-local cor- 
relations provided by AE^ 1 in vdW-DF. The figure also 
shows that the difference partly cancels, at short but 
not at intermediate separations, with the energy differ- 
ence between the exchange flavors of revPBE and PBE, 
SE X = AE™ vPBE - E PBE . Thus, for certain exchange 
functionals, adding an asymptotic atom-based account of 
non-local correlation can generate good binding values, 
yet our vdW-DF results indicate that this framework is 
not optimal for describing interactions at intermediate 
separations. 

Our results suggest that APPs and DFT-Ds could be 
improved at intermediate separations. A possible strat- 
egy is to replace the atomic separations r by an effective 
separation r— ro, where ro reflects the image-planes found 
for nanotubes and surfaces in Refs. [23ll7lT Keeping this 
(surface-physics) effect would increase the strength of the 
dispersion interactions at shorter separations. 

In summary, the results of this section suggests that 
an asymptotic atom-based pair potentials has a limited 
transferability over the full range of separations. Thus, 
for schemes using such a form, our results raises ques- 
tions on their ability to generate accurate results under 
broad condition (having multiple characteristic separa- 
tions), for instance involving phase transitions or pro- 
cesses that drive the system out off equilibrium, in pro- 
tein unfolding, in phase transitions, or simply for sys- 
tems which have competing interactions^ We note that 
there is no guarantee that vdW-DF, in its current form, 
can provide an accurate account under such broad con- 
ditions. A vdW-DF limitation is here exemplified by the 
likely underestimation of Cq coefficients for the platonic 
molecules. Nevertheless, we argue that non-local func- 
tionals, like vdW-DF, hold the most promise for deal- 
ing with molecular configurations under broad condi- 
tions. This is because an electron-based approach pro- 
vides a framework which naturally includes image-plane 
and multipole effects. It therefore holds the key to an ac- 
count which describe the variation in dispersive response 
over the full range of separations. 



VI. CONCLUSIONS 

For the three cage- like organic molecular crystals, hex- 
amine, cubane and dodecahedrane, vdW-DF gives lattice 
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parameters and cohesion energy that agree well with ex- 
periments, although for all three crystals the unit cell 
volumes are overestimated. A substantial sensitivity 
of lattice parameters and cohesion energy to the flavor 
of semi-local exchange signals that this overestimation 
stems mostly from the chosen form of exchange func- 
tional. 

We have also shown that, at short and intermediate 
separations, the full non-local correlation of vdW-DF is 
considerably larger than its corresponding atom-based 
asymptotic account. Notwithstanding that the asymp- 
totic account of vdW-DF likely needs improvement, this 
enhancement indicates that the asymptotic 1 /r e form of 
atomic pair potentials, by construction, can not give a 
transferable account over a large range of separations. 

This paper underlines the usefulness of studying sim- 
ple model sparse systems as molecular crystals to gain 
insight into methods intended for the study of sparse 
and supramolecular systems ~ Both DFT-D and vdW- 
DF benefit from such testing, because they are designed 
to be parameter-free, and an accurate account of molec- 



ular crystals would suggest an accurate account of more 
complex assemblies of similar molecules. The presented 
molecular crystals provide particularly accessible cases of 
extended systems and can therefore be used in conjunc- 
tion with future development of pair-potential methods 
and exchange-correlation functionals. 
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